##### Support Functions	

	### Compute Logit Gradient and Hessian

		# Reference: p.18 of https://cran.r-project.org/web/packages/mlogit/vignettes/mlogit.pdf
		
		function compute_nested_logit_gradient(betap,X,choices,probs,household_numbers,data_weights)
		
			lambda = last(betap);
			betap = betap[1:(length(betap)-1)];
			U = X * betap;
			
			inside_option_prob = 1 .- probs[uninsured_plans .== 1];
			inside_option_probs = inside_option_prob[household_numbers];
			
			ex_probs = probs ./ inside_option_probs;
			ex_probs[uninsured_plans .== 1] .= 0;
			
			probsdf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
			N_exchange_household = combine(groupby(probsdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			N_exchange = N_exchange_household[household_numbers];
			
			# Gradient wrt beta
			gradient_beta = (1/lambda) * (transpose(X) * (choices .* data_weights)) + 
							((lambda-1)/lambda) * (transpose(X) * (ex_probs .* I_ex_nest .* data_weights)) - 
							transpose(X) * (probs .* data_weights);
			
			# Gradient wrt lambda (This should be a scalar!)
			gradient_lambda = -(1/lambda^2) * (transpose(U) * (choices .* data_weights)) + 
							transpose(choices .* data_weights .* (1 .- uninsured_plans)) * log.(N_exchange) + 
							((1 - lambda)/lambda^2) * (transpose(U) * (ex_probs .* I_ex_nest .* data_weights)) .- 
							transpose(choices .* data_weights .* inside_option_probs) * log.(N_exchange) + 
							(1/lambda) * (transpose(U) * (probs .* data_weights));
			
			return(cat(dims=1,gradient_beta,gradient_lambda))
		end
		
		function compute_nested_logit_hessian(betap,X,choices,probs,household_numbers,data_weights)
		
			min_beta = betap;
			lambda = last(betap);
			betap = betap[1:(length(betap)-1)];
			U = X * betap;
			
			inside_option_prob = 1 .- probs[uninsured_plans .== 1];
			inside_option_probs = inside_option_prob[household_numbers];
			
			ex_probs = probs ./ inside_option_probs;
			ex_probs[uninsured_plans .== 1] .= 0;
			
			probsdf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								weighted_utility = probs[uninsured_plans .== 0] .* U[uninsured_plans .== 0],
								ex_weighted_utility = ex_probs[uninsured_plans .== 0] .* U[uninsured_plans .== 0];
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
			
			N_exchange_household = combine(groupby(probsdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			N_exchange = N_exchange_household[household_numbers];
		
			utility_sum_uni_household = combine(groupby(probsdf,:household_number), :weighted_utility => sum)[!,:weighted_utility_sum];
			utility_sum_uni = utility_sum_uni_household[household_numbers];
			utility_sum = copy(utility_sum_uni);
			utility_sum[uninsured_plans .== 1] .= 0;
		
			ex_utility_sum_uni_household = combine(groupby(probsdf,:household_number), :ex_weighted_utility => sum)[!,:ex_weighted_utility_sum];
			ex_utility_sum_uni = ex_utility_sum_uni_household[household_numbers];
			ex_utility_sum = copy(ex_utility_sum_uni);
			ex_utility_sum[uninsured_plans .== 1] .= 0;
			
			probX = DataFrame(probs .* X,:auto);
			probX[!,:household_number] = household_numbers;
			inner_sum_households = zeros(length(N_exchange_household),size(X)[2]);
			for k in 1:size(inner_sum_households)[2]
				inner_sum_households[:,k] = combine(groupby(probX,:household_number), k => sum)[!,2];
			end
			inner_sum = inner_sum_households[household_numbers,:];
			
			exprobX = DataFrame(ex_probs .* X,:auto);
			exprobX[!,:household_number] = household_numbers;
			ex_inner_sum_households = zeros(length(N_exchange_household),size(X)[2]);
			for k in 1:size(ex_inner_sum_households)[2]
				ex_inner_sum_households[:,k] = combine(groupby(exprobX,:household_number), k => sum)[!,2];
			end
			ex_inner_sum = ex_inner_sum_households[household_numbers,:];
			
			probX = nothing;
			exprobX = nothing;
			probsdf = nothing;
			N_exchange_household = nothing;
			utility_sum_uni_household = nothing;
			ex_utility_sum_uni_household = nothing;
			inner_sum_households = nothing;
			ex_inner_sum_households = nothing;
			GC.gc()
			
			# Partial derivatives we need to build Hessian
			partial_ex_prob_partial_beta = (1/lambda) .* (ex_probs .* I_ex_nest) .* (X .- ex_inner_sum);
			partial_prob_partial_beta = (1/lambda) .* (X .* probs) .+ (lambda - 1)/lambda .* (probs .* ex_inner_sum) .- (probs .* inner_sum);
			partial_logNexchange_partial_beta = (1/lambda) .* ex_inner_sum;
			partial_inside_prob_partial_beta = (1 .- inside_option_probs) .* inner_sum;
			
			
			#partial_prob_partial_beta = probs .* X;
			#partial_prob_partial_beta .*= 1/lambda;
			#partial_prob_partial_beta += (lambda - 1)/lambda * (probs .* ex_inner_sum);
			ex_inner_sum = nothing;
			GC.gc()
			#partial_prob_partial_beta -= probs .* inner_sum;
			inner_sum = nothing;
			GC.gc()
			
			partial_logNexchange_partial_lambda = -(1/lambda^2) * ex_utility_sum;
			partial_logNexchange_partial_lambda_uni = -(1/lambda^2) * ex_utility_sum_uni;
			partial_ex_prob_partial_lambda = -(1/lambda^2) * ex_probs .* (U - ex_utility_sum);
			partial_inside_prob_partial_lambda = -(1/lambda) * (1 .- inside_option_probs) .* (utility_sum) + 
				(1 .- inside_option_probs) .* inside_option_probs .* log.(N_exchange) ;
			partial_inside_prob_partial_lambda_uni = -(1/lambda) * (1 .- inside_option_probs) .* (utility_sum_uni) + 
				(1 .- inside_option_probs) .* inside_option_probs .* log.(N_exchange) ;
			partial_prob_partial_lambda = -(1/lambda^2) * (U .* probs) - ((lambda-1)/lambda^2) * probs .* ex_utility_sum + 
				(probs .* log.(N_exchange)) - probs .* (-1/lambda * utility_sum + inside_option_probs .* log.(N_exchange));
			
			# Now we build the Hessian
			
			hess = zeros(length(min_beta),length(min_beta));	
				
				# Covariate part of hessian
				hess[1:length(betap),1:length(betap)] = ((lambda - 1)/(lambda)) .*
					transpose(data_weights .* X) * (partial_ex_prob_partial_beta .* I_ex_nest) .-
					transpose(data_weights .* X) * partial_prob_partial_beta;	
				
				# Last row and column of hessian (except right corner)
				hess[1:length(betap),size(hess)[1]] = -(1/lambda^2) .* (transpose(data_weights .* X) * choices) .+
					transpose(partial_logNexchange_partial_beta) * (data_weights .* choices .* (1 .- uninsured_plans)) .+ 
					((1 - lambda)/lambda^2) .* (transpose(partial_ex_prob_partial_beta .* I_ex_nest) * (data_weights .* U) .+ 
						transpose(data_weights .* X) * (ex_probs .* I_ex_nest)) .-
					transpose(partial_inside_prob_partial_beta .* log.(N_exchange) .+
						inside_option_probs .* partial_logNexchange_partial_beta) * (data_weights .* choices) .+ 
					(1/lambda) .* (transpose(partial_prob_partial_beta) * (data_weights .*  U) .+ 
						transpose(data_weights .* X) * probs) ;					
									
				hess[size(hess)[1],1:length(betap)] = hess[1:length(betap),size(hess)[1]];					
									
				# Far right corner of Hessian
				hess[size(hess)[1],size(hess)[1]] = ((2/lambda^3) * (transpose(U) * (data_weights .* choices)) +
								sum((data_weights .* choices .* (1 .- uninsured_plans)) .* partial_logNexchange_partial_lambda) +
								((1 - lambda)/lambda^2) * (transpose(U) * (data_weights .* partial_ex_prob_partial_lambda .* I_ex_nest)) + 
									(lambda-2)/lambda^3 * (transpose(U) * (data_weights .* ex_probs .* I_ex_nest)) - 
								(transpose(data_weights .* choices .* partial_inside_prob_partial_lambda_uni) * log.(N_exchange) + 
									transpose(data_weights .* choices .* inside_option_probs) * partial_logNexchange_partial_lambda_uni) +
								(1/lambda) * (transpose(U) * (data_weights .* partial_prob_partial_lambda)) - 
									(1/lambda^2) * (transpose(U) * (data_weights .* probs)))[1];
			
			return(hess)
		end
		
		function compute_logit_gradient(X,choices,probs,data_weights)
			return(transpose(X) * ((choices .- probs).* data_weights))
		end
		
		function compute_logit_hessian(betap,X,choices,probs,household_numbers,data_weights)
		
			probX = DataFrame(probs .* X,:auto);
			probX[!,:household_number] = household_numbers;
			inner_sum_households = zeros(size(households)[1],size(X)[2]);
			for k in 1:size(inner_sum_households)[2]
				inner_sum_households[:,k] = combine(groupby(probX,:household_number), k => sum)[!,2];
			end
			inner_sum = inner_sum_households[household_numbers,:];
			
			partial_prob_partial_beta = (X .- inner_sum) .* probs;
			hess = -transpose(data_weights .* X) * partial_prob_partial_beta;	
				
			return(hess)
		end
		
		function compute_rp_logit_gradient(X,choices,probabilities,prob_matrix,data_weights)
		
			Xupdated = copy(X);
			gradient = zeros(length(all_covariates));
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				probsdf = DataFrame(household_number = household_numbers);
				probX = Xupdated .* probsn; # IJ * K - we need to sum within household for each covariate
				Xprob = zeros(length(probabilities),length(all_covariates));
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = probX[:,k];
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					Xprob[:,k] = probsum[household_numbers];
				end
				gradient += 1/ns .* (transpose(Xupdated .- Xprob) * (data_weights .* choices .* probsn ./ probabilities));
			end
		
			return(gradient)
		end
		
		function compute_rp_logit_hessian(X,choices,probabilities,prob_matrix,data_weights)
		
			hess = zeros(length(all_covariates),length(all_covariates)); # K * K
		
			# Get average probability
			
			probsdf = DataFrame(zeros(length(household_numbers),length(all_covariates)),Symbol.(all_covariates));
			probsdf[!,:household_number] = household_numbers;
			partial_averageprob_partial_beta_matrix = zeros(length(choices),length(all_covariates));
			for n in 1:ns
				Xupdated = copy(X);
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] .*= V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				
				probsdf[!,Symbol.(all_covariates)] = Xupdated .* probsn;
				probsum = Matrix(combine(groupby(probsdf,:household_number), Symbol.(all_covariates) .=> sum; renamecols=false)[household_numbers,Symbol.(all_covariates)]);
				partial_averageprob_partial_beta_matrix += 1/ns * ((Xupdated .- probsum) .* probsn);
			end
		
			for n in 1:ns
				Xupdated = copy(X);
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] .*= V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				
				probsdf[!,Symbol.(all_covariates)] = Xupdated .* probsn;
				probsum = Matrix(combine(groupby(probsdf,:household_number), Symbol.(all_covariates) .=> sum; renamecols=false)[household_numbers,Symbol.(all_covariates)]);
				partial_probn_partial_beta = probsn .* (Xupdated .- probsum); 
				partial_probratio_partial_beta = (partial_probn_partial_beta .* probabilities  .- partial_averageprob_partial_beta_matrix .* probsn)./ (probabilities).^2;
					
				for kp in 1:length(all_covariates)
					probsdf[!,Symbol.(all_covariates)] = Xupdated .* partial_probn_partial_beta[:,kp];
					partial_probsum = Matrix(combine(groupby(probsdf, :household_number), Symbol.(all_covariates) .=> sum; renamecols=false)[household_numbers,Symbol.(all_covariates)]);
					hess[:,kp] += 1/ns .* (transpose(Xupdated .- probsum) * (data_weights .* choices .* partial_probratio_partial_beta[:,kp]) .- 
						transpose(partial_probsum) * (data_weights .* choices .* probsn ./ probabilities));
				end
			end
		
			return(hess)
		end		
			
		
		function compute_rp_logit_hessian_old(X,choices,probabilities,prob_matrix,data_weights)
		
			hess = zeros(length(all_covariates),length(all_covariates)); # K * K
		
			# Get probability inner sum and average probability
			# It would be less computationally intensive to save output in 3D array, but not enough memory to do this!
		
			probsdf = DataFrame(household_number = household_numbers);
			partial_averageprob_partial_beta_matrix = zeros(length(choices),length(all_covariates));
			for n in 1:ns
				Xupdated = copy(X);
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] .*= V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					partial_averageprob_partial_beta_matrix[:,k] += 1/ns * (probsn .* (Xupdated[:,k] .- probsum[household_numbers]));
				end
			end
		
			for n in 1:ns
				Xupdated = copy(X);
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] .*= V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				for k in 1:length(all_covariates)
					Xk = Xupdated[:,k];
					probsdf[!,:probXk] = Xk .* probsn;
					probsumk = combine(groupby(probsdf,:household_number), :probXk => sum)[!,:probXk_sum];
					for kp in k:length(all_covariates)
						Xkp = Xupdated[:,kp];
						probsdf[!,:probXkp] = Xkp .* probsn;
						probsumkp = combine(groupby(probsdf,:household_number), :probXkp => sum)[!,:probXkp_sum];
						partial_probn_partial_beta = probsn .* (Xkp .- probsumkp[household_numbers]);
				
						partial_probratio_partial_beta = 
							(partial_probn_partial_beta .* probabilities  .- partial_averageprob_partial_beta_matrix[:,kp] .* probsn)./ (probabilities).^2;
				
						probsdf[!,:partialprobX] = Xk .* partial_probn_partial_beta;
						partial_probsum = combine(groupby(probsdf,:household_number), :partialprobX => sum)[!,:partialprobX_sum];
						partialXprob = partial_probsum[household_numbers]; 
				
						hess[k,kp] += 1/ns .* (sum(data_weights .* choices .* (Xk .- probsumk[household_numbers]) .* partial_probratio_partial_beta) .-
							sum(partialXprob .* data_weights .* choices .* probsn ./ probabilities)); 
						hess[kp,k] = hess[k,kp];	
					end
				end
			end
		
			return(hess)
		end
		
		
		
		function compute_rp_nested_logit_gradient(betap,X,choices,probabilities,prob_matrix,data_weights,household_numbers)
		
			Xupdated = copy(X);
			gradient = zeros(length(all_covariates)+1);
		
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
		
			for n in 1:ns
				
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,n]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				# Gradient wrt beta
				probsdf = DataFrame(household_number = household_numbers);
				gradient_beta = zeros(length(all_covariates));
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,k] .* ex_probs;
					
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					 
					gradient_beta[k] =  sum(((1/lambda) .* Xupdated[:,k] .+ ((lambda-1)/lambda) .* Xexprobk .* I_ex_nest .- Xprobk) .* choices .* data_weights .* probsn ./ probabilities);
				end
				
				# Gradient wrt lambda (This should be a scalar!)
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
				Uprob = probUsum[household_numbers];
				Uexprob = exprobUsum[household_numbers];		
				
				probsdf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
				N_exchange_household = combine(groupby(probsdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				N_exchange = N_exchange_household[household_numbers];
				
				gradient_lambda = sum((-(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
								((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .-
								inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob) .*
								choices .* data_weights .* probsn ./ probabilities);

				gradient += 1/ns .* cat(dims=1,gradient_beta,gradient_lambda);
			end
			
			return(gradient)
		end
		
		function compute_rp_nested_logit_hessian(betap,X,choices,probabilities,prob_matrix,data_weights,household_numbers)
		
			hess = zeros(length(all_covariates)+1,length(all_covariates)+1); # K * K
			Xupdated = copy(X);
			
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
			
			probsdf = DataFrame(household_number = household_numbers);
			partial_averageprob_partial_beta_matrix = zeros(length(choices),length(all_covariates));
			partial_averageprob_partial_lambda = zeros(length(choices));
			
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,n]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
				Uprob = probUsum[household_numbers];
				Uexprob = exprobUsum[household_numbers];	
				
				Ndf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
				N_exchange_household = combine(groupby(Ndf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				N_exchange = N_exchange_household[household_numbers];
				
				nested_lambda_multiplier = -(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
					((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .- inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob;
				
				partial_probn_partial_lambda = probsn .* nested_lambda_multiplier;
				partial_averageprob_partial_lambda += 1/ns .* partial_probn_partial_lambda;	
				
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,k] .* ex_probs;
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					
					partial_probn_partial_beta = probsn .* ((1/lambda) .* Xupdated[:,k] .+ ((lambda-1)/lambda) .* Xexprobk .* I_ex_nest .- Xprobk);
					partial_averageprob_partial_beta_matrix[:,k] += 1/ns .* partial_probn_partial_beta;
				end			
			end
			
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,n]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
				Uprob = probUsum[household_numbers];
				Uexprob = exprobUsum[household_numbers];	
				
				Ndf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
				N_exchange_household = combine(groupby(Ndf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				N_exchange = N_exchange_household[household_numbers];
					
				nested_lambda_multiplier = -(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
					((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .- inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob;
				partial_probn_partial_lambda = probsn .* nested_lambda_multiplier;
				partial_ex_prob_partial_lambda = -(1/lambda^2) .* ex_probs .* (U .- Uexprob);
				partial_logNexchange_partial_lambda = -(1/lambda^2) .* Uexprob;
				partial_inside_prob_partial_lambda = -(1/lambda) .* (1 .- inside_option_probs) .* Uprob .+ (1 .- inside_option_probs) .* inside_option_probs .* log.(N_exchange);
				
				probsdf[!,:partialprobUlambda] = U .* partial_probn_partial_lambda;
				probsdf[!,:partialexprobUlambda] = U .* partial_ex_prob_partial_lambda;
				partial_probsumUlambda = combine(groupby(probsdf,:household_number), :partialprobUlambda => sum)[!,:partialprobUlambda_sum];
				partial_exprobsumUlambda = combine(groupby(probsdf,:household_number), :partialexprobUlambda => sum)[!,:partialexprobUlambda_sum];
				partial_Uprob_partial_lambda = partial_probsumUlambda[household_numbers]; 
				partial_Uexprob_partial_lambda = partial_exprobsumUlambda[household_numbers]; 			
			
				partial_probratio_partial_lambda = (partial_probn_partial_lambda .* probabilities  .- partial_averageprob_partial_lambda .* probsn)./ (probabilities).^2;
				partial_nested_lambda_multiplier_partial_lambda = (2/lambda^3) .* U .+ 
					(1 .- uninsured_plans) .* partial_logNexchange_partial_lambda .+ 
					((1 - lambda)/lambda^2) .* partial_Uexprob_partial_lambda .* I_ex_nest .+ (lambda-2)/lambda^3 .* Uexprob .* I_ex_nest .-
					(partial_inside_prob_partial_lambda .* log.(N_exchange) .+ inside_option_probs .* partial_logNexchange_partial_lambda) .+ 
					(1/lambda) .* partial_Uprob_partial_lambda .* I_ex_nest .- (1/lambda^2) .* Uprob;
				
				# Far right corner of Hessian
				hess[size(hess)[1],size(hess)[1]] += 1/ns .* sum((nested_lambda_multiplier .* partial_probratio_partial_lambda .+
							partial_nested_lambda_multiplier_partial_lambda .* probsn ./ probabilities) .* choices .* data_weights);
				
				for k in 1:length(all_covariates)
					probsdf[!,:probXk] = Xupdated[:,k] .* probsn;
					probsdf[!,:exprobXk] = Xupdated[:,k] .* ex_probs;
					probsumk = combine(groupby(probsdf,:household_number), :probXk => sum)[!,:probXk_sum];
					exprobsumk = combine(groupby(probsdf,:household_number), :exprobXk => sum)[!,:exprobXk_sum];
					nested_multiplierk = ((1/lambda) .* Xupdated[:,k] .+ ((lambda-1)/lambda) .* exprobsumk[household_numbers] .* I_ex_nest .- probsumk[household_numbers]);
					
					partial_probn_partial_betak = probsn .* nested_multiplierk;
					partial_ex_prob_partial_betak = (1/lambda) .* ex_probs .* (Xupdated[:,k] .- exprobsumk[household_numbers]);
					partial_inside_prob_partial_betak = (1 .- inside_option_probs) .* probsumk[household_numbers];
					partial_logNexchange_partial_betak = (1/lambda) .* exprobsumk[household_numbers];
						
					probsdf[!,:partialprobU] = U .* partial_probn_partial_betak;
					probsdf[!,:partialexprobU] = U .* partial_ex_prob_partial_betak;
					partial_probsumU = combine(groupby(probsdf,:household_number), :partialprobU => sum)[!,:partialprobU_sum];
					partial_exprobsumU = combine(groupby(probsdf,:household_number), :partialexprobU => sum)[!,:partialexprobU_sum];
					partialUprob = partial_probsumU[household_numbers]; 
					partialUexprob = partial_exprobsumU[household_numbers]; 
					partial_Uprob_partial_betak = partialUprob .* I_ex_nest .+ probsumk[household_numbers];
					partial_Uexprob_partial_betak = (partialUexprob .+ exprobsumk[household_numbers]) .* I_ex_nest;
							
					partial_probratio_partial_betak = (partial_probn_partial_betak .* probabilities  .- partial_averageprob_partial_beta_matrix[:,k] .* probsn)./ (probabilities).^2;
					partial_nested_lambda_multiplier_partial_betak = -(1/lambda^2) .* Xupdated[:,k] .+ 
						(1 .- uninsured_plans) .* partial_logNexchange_partial_betak .+ 
						((1 - lambda)/lambda^2) .* partial_Uexprob_partial_betak .- 
						(partial_inside_prob_partial_betak .* log.(N_exchange) .+ inside_option_probs .* partial_logNexchange_partial_betak) .+ 
						(1/lambda) .* partial_Uprob_partial_betak;
					
					# Bottom row and last column (not far right corner)
					hess[size(hess)[1],k] += 1/ns .* sum((nested_lambda_multiplier .* partial_probratio_partial_betak .+
						partial_nested_lambda_multiplier_partial_betak .* probsn ./ probabilities) .* choices .* data_weights);	
						
					for kp in k:length(all_covariates)
						probsdf[!,:probXkp] = Xupdated[:,kp] .* probsn;
						probsdf[!,:exprobXkp] = Xupdated[:,kp] .* ex_probs;
						probsumkp = combine(groupby(probsdf,:household_number), :probXkp => sum)[!,:probXkp_sum];
						exprobsumkp = combine(groupby(probsdf,:household_number), :exprobXkp => sum)[!,:exprobXkp_sum];
						nested_multiplierkp = (1/lambda) .* Xupdated[:,kp] .+ ((lambda-1)/lambda) .* exprobsumkp[household_numbers] .* I_ex_nest .- probsumkp[household_numbers];
						
						partial_probn_partial_betakp = probsn .* nested_multiplierkp;
						partial_ex_prob_partial_betakp = (1/lambda) .* ex_probs .* (Xupdated[:,kp] .- exprobsumkp[household_numbers]);
						
						probsdf[!,:partialprobX] = Xupdated[:,k] .* partial_probn_partial_betakp;
						probsdf[!,:partialexprobX] = Xupdated[:,k] .* partial_ex_prob_partial_betakp;
						partial_probsum = combine(groupby(probsdf,:household_number), :partialprobX => sum)[!,:partialprobX_sum];
						partial_exprobsum = combine(groupby(probsdf,:household_number), :partialexprobX => sum)[!,:partialexprobX_sum];
						partialXprob = partial_probsum[household_numbers]; 
						partialXexprob = partial_exprobsum[household_numbers]; 
						
						partial_probratio_partial_betakp = (partial_probn_partial_betakp .* probabilities  .- partial_averageprob_partial_beta_matrix[:,kp] .* probsn)./ (probabilities).^2;
						partial_nested_multiplier_partial_betakp = ((lambda - 1)/lambda .* partialXexprob .- partialXprob) .* I_ex_nest;
							
						# Covariate portion of Hessian
						hess[k,kp] += 1/ns .* sum((nested_multiplierk .* partial_probratio_partial_betakp .+
							partial_nested_multiplier_partial_betakp .* probsn ./ probabilities) .* choices .* data_weights);
						hess[kp,k] = hess[k,kp];
					end
				end
			end
								
			hess[1:length(betap),size(hess)[1]] = vec(hess[size(hess)[1],1:length(betap)]);					
			
			return(hess)
		end
		
		function compute_rp_logit_seq_gradient(X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights)
		
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			
			probsdf = DataFrame(household_number = household_numbers);
			gradient = zeros(length(all_covariates));
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				sequence_probsn = vec(sequence_prob_matrix[:,n]);
				probX = Xupdated .* probsn;
				
				Xseq = zeros(length(hseq_weights),length(all_covariates));
				Xseq_prob = zeros(length(sequence_probs),length(all_covariates));
				
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = probX[:,k];
					hseq_probs_updated[!,:X] = Xupdated[observed_choice_indices,k]; 
					hseq_probs_updated[!,:Xseq_prob] = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					
					Xseq[:,k] = combine(groupby(hseq_probs_updated,:household_id), :X => sum)[!,:X_sum];
					Xseq_prob[:,k] = combine(groupby(hseq_probs_updated,:household_id), :Xseq_prob => sum)[!,:Xseq_prob_sum];
				end
				gradient += 1/ns .* (transpose(Xseq .- Xseq_prob) * (hseq_weights .* sequence_probsn ./ sequence_probs));
			end
		
			return(gradient)
		end
		
		function compute_rp_nested_logit_seq_gradient(betap,X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights)
		
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			gradient = zeros(length(all_covariates)+1);
		
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
		
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,n]);
				sequence_probsn = vec(sequence_prob_matrix[:,n]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				hseq_probs_updated[!,:inside_option_prob] = inside_option_prob; 
				inside_option_seq = combine(groupby(hseq_probs_updated,:household_id), :inside_option_prob => prod)[!,:inside_option_prob_prod];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				# Gradient wrt beta
				probsdf = DataFrame(household_number = household_numbers);
				gradient_beta = zeros(length(all_covariates));
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,k] .* ex_probs;
					
					hseq_probs_updated[!,:Xseq] = Xupdated[observed_choice_indices,k]; 
					hseq_probs_updated[!,:Xseq_prob] = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					hseq_probs_updated[!,:Xseq_exprob] = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum] .* I_ex_nest_household;
					
					Xseqk = combine(groupby(hseq_probs_updated,:household_id), :Xseq => sum)[!,:Xseq_sum];
					Xseq_probk = combine(groupby(hseq_probs_updated,:household_id), :Xseq_prob => sum)[!,:Xseq_prob_sum];
					Xseq_exprobk = combine(groupby(hseq_probs_updated,:household_id), :Xseq_exprob => sum)[!,:Xseq_exprob_sum];
					
					nested_multiplier = (1/lambda) .* Xseqk .+ ((lambda-1)/lambda) .* Xseq_exprobk .- Xseq_probk;
					
					gradient_beta[k] =  sum(nested_multiplier .* hseq_weights .* sequence_probsn ./ sequence_probs);
				end
				
				# Gradient wrt lambda (This should be a scalar!)
					# NOTE: I do not get the right answer.  May consider using built-in gradient function instead for lamda partial
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probsdf[!,:N_exchange] = (1 .- uninsured_plans) .* exp.(U ./ lambda);
				hseq_probs_updated[!,:Useq] = U[observed_choice_indices];
				hseq_probs_updated[!,:Useq_prob] = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				hseq_probs_updated[!,:Useq_exprob] = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum] .* I_ex_nest_household;
				hseq_probs_updated[!,:Nseq_exchange] = combine(groupby(probsdf,:household_number), :N_exchange => sum)[!,:N_exchange_sum];
				#hseq_probs_updated[!,:Nseq_exchange] = exp.(U[observed_choice_indices] ./ lambda);
				#hseq_probs_updated[!,:logNseq_exchange] .= 0.0;
				#update_indices = findall(hseq_probs_updated[!,:Nseq_exchange] .> 0);
				hseq_probs_updated[!,:logNseq_exchange] = log.(hseq_probs_updated[!,:Nseq_exchange]);
				
				Useq = combine(groupby(hseq_probs_updated,:household_id), :Useq => sum)[!,:Useq_sum];
				Useq_prob = combine(groupby(hseq_probs_updated,:household_id), :Useq_prob => sum)[!,:Useq_prob_sum];
				Useq_exprob = combine(groupby(hseq_probs_updated,:household_id), :Useq_exprob => sum)[!,:Useq_exprob_sum];
				logN_exchange_seq = combine(groupby(hseq_probs_updated,:household_id), :logNseq_exchange => sum)[!,:logNseq_exchange_sum];
								
				nested_lambda_multiplier = -(1/lambda^2) .* Useq .+ logN_exchange_seq .+  
					((1 - lambda)/lambda^2) .* Useq_exprob .- inside_option_seq .* logN_exchange_seq .+ (1/lambda) .* Useq_prob; 
				gradient_lambda = sum(nested_lambda_multiplier .* hseq_weights .* sequence_probsn ./ sequence_probs);

				gradient_lambda = derivative(calculate_partial_nested_prob_seq_partial_lambda,lambda);

				gradient += 1/ns .* cat(dims=1,gradient_beta,gradient_lambda);
			end
			
			return(gradient)
		end
		
		function compute_rp_logit_seq_hessian(X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights)
		
			hess = zeros(length(all_covariates),length(all_covariates)); # K * K
			
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			
			probsdf = DataFrame(household_number = household_numbers);
			partial_averageseqprob_partial_beta_matrix = zeros(length(hseq_weights),length(all_covariates));
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				sequence_probsn = vec(sequence_prob_matrix[:,n]);
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					hseq_probs_updated[!,:X] = Xupdated[observed_choice_indices,k]; 
					hseq_probs_updated[!,:Xseq_prob] = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					Xseqk = combine(groupby(hseq_probs_updated,:household_id), :X => sum)[!,:X_sum];
					Xseqk_prob = combine(groupby(hseq_probs_updated,:household_id), :Xseq_prob => sum)[!,:Xseq_prob_sum];
					partial_averageseqprob_partial_beta_matrix[:,k] += 1/ns .* (sequence_probsn .* (Xseqk .- Xseqk_prob));
				end
			end
		
			for n in 1:ns
				Xupdated = copy(X);
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				sequence_probsn = vec(sequence_prob_matrix[:,n]);
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					hseq_probs_updated[!,:X] = Xupdated[observed_choice_indices,k]; 
					hseq_probs_updated[!,:Xseq_prob] = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					Xseqk = combine(groupby(hseq_probs_updated,:household_id), :X => sum)[!,:X_sum];
					Xseqk_prob = combine(groupby(hseq_probs_updated,:household_id), :Xseq_prob => sum)[!,:Xseq_prob_sum];
					
					for kp in k:length(all_covariates)
						probsdf[!,:probXkp] = Xupdated[:,kp] .* probsn;
						hseq_probs_updated[!,:Xkp] = Xupdated[observed_choice_indices,kp]; 
						hseq_probs_updated[!,:Xseqkp_prob] = combine(groupby(probsdf,:household_number), :probXkp => sum)[!,:probXkp_sum];
						Xseqkp = combine(groupby(hseq_probs_updated,:household_id), :Xkp => sum)[!,:Xkp_sum];
						Xseqkp_prob = combine(groupby(hseq_probs_updated,:household_id), :Xseqkp_prob => sum)[!,:Xseqkp_prob_sum];
					
						partial_probn_partial_beta = probsn .* (Xupdated[:,kp] .- hseq_probs_updated[household_numbers,:Xseqkp_prob]);
						partial_seqprobn_partial_betakp = sequence_probsn .* (Xseqkp .- Xseqkp_prob)
						
						partial_probratio_partial_betakp = (sequence_probs .* partial_seqprobn_partial_betakp  .- 
							sequence_probsn .* partial_averageseqprob_partial_beta_matrix[:,kp]) ./ (sequence_probs).^2;
					
						probsdf[!,:partialprobseqX] = Xupdated[:,k] .* partial_probn_partial_beta;
						hseq_probs_updated[!,:partialprobseqX] = combine(groupby(probsdf,:household_number), :partialprobseqX => sum)[!,:partialprobseqX_sum];
						Xseqkp_prob_partial = combine(groupby(hseq_probs_updated,:household_id), :partialprobseqX => sum)[!,:partialprobseqX_sum];
					
						hess[k,kp] += 1/ns * (sum(hseq_weights .* (Xseqk .- Xseqk_prob) .* partial_probratio_partial_betakp) .-
							sum(hseq_weights .* sequence_probsn ./ sequence_probs .* Xseqkp_prob_partial)); 
						hess[kp,k] = hess[k,kp];
					end
				end
			end
			
			return(hess)
		end
		
		
	### Simulate Probabilities
	
		function simulate_nested_probabilities(betap,Xins,household_numbers_ins,uninsured_plans) 
			lambda = last(betap);
			betap = betap[1:(length(betap)-1)];

			expdf = DataFrame(household_number = household_numbers_ins,expvarall = exp.((Xins * betap)/lambda));
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			
			probs = zeros(length(uninsured_plans));
			probs[uninsured_plans .== 0] = expdf[!,:expvarall] .* (expsum[household_numbers_ins]).^(lambda-1) ./ (1 .+ (expsum[household_numbers_ins]).^(lambda));
			probs[uninsured_plans .== 1] = 1  ./ (1 .+ (expsum).^(lambda));
			
			return(probs)
		end	
		
		function simulate_logit_probabilities(betap,Xins,household_numbers_ins,uninsured_plans) 
			expdf = DataFrame(household_number = household_numbers_ins,expvarall = exp.(Xins * betap));
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			
			probs = zeros(length(uninsured_plans));
			probs[uninsured_plans .== 0] = expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers_ins]);
			probs[uninsured_plans .== 1] = 1  ./ (1 .+ expsum);
			
			return(probs)
		end	

		function simulate_logit_probabilities_smooth(betap,Xins,household_numbers_ins,uninsured_plans,tau) 
			expdf = DataFrame(household_number = household_numbers_ins,expvarall = exp.((Xins * betap) ./ tau));
			expsum = combine(groupby(expdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			
			probs = zeros(length(uninsured_plans));
			probs[uninsured_plans .== 0] = expdf[!,:expvarall] ./ (1 .+ expsum[household_numbers_ins]);
			probs[uninsured_plans .== 1] = 1  ./ (1 .+ expsum);
			
			return(probs)
		end	
	
	### Calculate Search Probabilities

		function calculate_search_probabilities(beta_search,Xsearch)
			household_search_probs = 1 ./ (1 .+ exp.(Xsearch * beta_search));
			return(household_search_probs[household_numbers])
		end
		
	### Calculate Log Likelihood
	
		function calculate_log_likelihood(choices,probs,data_weights)
			return((transpose(choices .* data_weights) * log.(probs))[1])
		end
		
		function calculate_sequence_log_likelihood(sequence_probs,hseq_weights)
			return(transpose(hseq_weights) * log.(sequence_probs))
		end
		
		function calculate_inattention_log_likelihood(choices,probs,search_probs)
			return((transpose(choices .* entrant_weights) * log.(probs))[1] .+ 
					(transpose(choices .* switcher_weights) * log.(search_probs .* probs))[1] .+ 
					(transpose(choices .* nonswitcher_weights) * log.((1 .- search_probs) .+ search_probs .* probs))[1])
		end
	
		function calculate_logit_log_likelihood_for_hessian(min_beta)
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			return((transpose(choices .* data_weights) * log.(probabilities))[1])
		end
	
		function calculate_sequence_log_likelihood_for_test(min_beta)
		
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				hseq_probs_updated[!,:prob] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],
					household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
				sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod];
			end
			sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
		
			return(transpose(hseq_weights) * log.(sequence_probs))
		end
	
		function calculate_nested_logit_log_likelihood_for_hessian(min_beta)
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			return((transpose(choices .* data_weights) * log.(probabilities))[1])
		end
	
		function calculate_rplogit_probability_gradient(min_beta)
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			return(probabilities[1])
		end
	
		function calculate_rp_nested_logit_probability_gradient(min_beta)
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			return(probabilities[1])
		end
	
		function calculate_rp_nested_logit_insideprob_gradient(min_beta)
			Xupdated = copy(X);
			Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
			
			probsn = simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
			inside_option_probs = inside_option_prob[household_numbers];
			
			return(inside_option_probs[100])
		end
	
		function calculate_rp_nested_logit_exprob_gradient(min_beta)
			Xupdated = copy(X);
			Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
			
			probsn = simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
			inside_option_probs = inside_option_prob[household_numbers];
			ex_probs = probsn ./ inside_option_probs;
			ex_probs[uninsured_plans .== 1] .= 0;
			
			return(ex_probs[100])
		end
		
		function calculate_rp_nested_logit_logNexchange_gradient(betap)
			Xupdated = copy(X);
			Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
			
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
			U = Xupdated * betanew;
		
			Ndf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
							expvarall = exp.(U[uninsured_plans .== 0]/lambda));
			N_exchange_household = combine(groupby(Ndf,:household_number), :expvarall => sum)[!,:expvarall_sum];
			N_exchange = N_exchange_household[household_numbers];
			logN_exchange = log.(N_exchange);
		
			return(logN_exchange[1])
		end
		
		function calculate_rp_probratio_gradient(betap)
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
			probsn = vec(prob_matrix[:,1]);
			probratio = probsn ./ probabilities;
			
			return(probratio[1])
		end
		
		function calculate_rp_nested_multiplier_gradient(betap)
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			Xupdated = copy(X);
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
		
			#for n in 1:ns
				
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,1];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,1]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				# Gradient wrt beta
				probsdf = DataFrame(household_number = household_numbers);
				#for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,1] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,1] .* ex_probs;
					
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					 
					nested_multiplier = (1/lambda) .* Xupdated[:,1] .+ ((lambda-1)/lambda) .* Xexprobk .* I_ex_nest .- Xprobk;
				#end
				
			#end
			
			return(nested_multiplier[1000000])
		end
		
		function calculate_rp_nested_lambda_multiplier_gradient(betap)
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			Xupdated = copy(X);
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
		
			#for n in 1:ns
				
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,1];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,1]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				probsdf = DataFrame(household_number = household_numbers);
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
				Uprob = probUsum[household_numbers];
				Uexprob = exprobUsum[household_numbers];	
				
				Ndf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
				N_exchange_household = combine(groupby(Ndf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				N_exchange = N_exchange_household[household_numbers];
					
				# Gradient wrt beta
				
				#for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,1] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,1] .* ex_probs;
					
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					 
					nested_lambda_multiplier = -(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
						((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .- inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob;
				#end
				
			#end
			
			return(nested_lambda_multiplier[1000000])
		end
		
		
		
		
		function calculate_rp_logit_gradient_for_test(min_beta)
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			gradient = 0.0;
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				probsn = vec(prob_matrix[:,n]);
				probsdf = DataFrame(household_number = household_numbers);
				probsdf[!,:probX] = Xupdated[:,1] .* probsn; # IJ * 1 - we need to sum within household for each covariate
				probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
				gradient += 1/ns .* (transpose(Xupdated[:,1] .- probsum[household_numbers]) * (data_weights .* choices .* probsn ./ probabilities));
			end
		
			return(gradient)
		end
	
		function calculate_rp_nested_logit_gradient_for_test(betap)
		
			beta_flag = false;
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			Xupdated = copy(X);
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
			
			#for n in 1:ns
				
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,1];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,1]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				# Gradient wrt beta
				probsdf = DataFrame(household_number = household_numbers);
				if beta_flag 
				#for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,1] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,1] .* ex_probs;
					
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					 
					nested_multiplier = (1/lambda) .* Xupdated[:,1] .+ ((lambda-1)/lambda) .* Xexprobk .* I_ex_nest .- Xprobk;
					 
					gradient = sum(nested_multiplier .* probsn ./ probabilities .* choices .* data_weights );
				#end
				else 
					# Gradient wrt lambda (This should be a scalar!)
					Ndf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
									expvarall = exp.(U[uninsured_plans .== 0]/lambda));
					N_exchange_household = combine(groupby(Ndf,:household_number), :expvarall => sum)[!,:expvarall_sum];
					N_exchange = N_exchange_household[household_numbers];
					
					probsdf[!,:probU] = U .* probsn;
					probsdf[!,:exprobU] = U .* ex_probs;
					probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
					exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
					Uprob = probUsum[household_numbers];
					Uexprob = exprobUsum[household_numbers];		
					
					probsdf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
									expvarall = exp.(U[uninsured_plans .== 0]/lambda));
					N_exchange_household = combine(groupby(probsdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
					N_exchange = N_exchange_household[household_numbers];
					
					gradient = sum((-(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
									((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .-
									inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob) .*
									choices .* data_weights .* probsn ./ probabilities);
				end
			#end
			
			return(gradient)
		end
	
		function calculate_rp_nested_logit_gradient!(grad,betap)
		
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			Xupdated = copy(X);
			gradient = zeros(length(all_covariates)+1);
			lambda = last(betap);
			betanew = betap[1:(length(betap)-1)];
		
			for n in 1:ns
				
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				U = Xupdated * betanew;
				
				probsn = vec(prob_matrix[:,n]);
				
				inside_option_prob = 1 .- probsn[uninsured_plans .== 1];
				inside_option_probs = inside_option_prob[household_numbers];
				
				ex_probs = probsn ./ inside_option_probs;
				ex_probs[uninsured_plans .== 1] .= 0;
				
				# Gradient wrt beta
				probsdf = DataFrame(household_number = household_numbers);
				gradient_beta = zeros(length(all_covariates));
				for k in 1:length(all_covariates)
					probsdf[!,:probX] = Xupdated[:,k] .* probsn;
					probsdf[!,:exprobX] = Xupdated[:,k] .* ex_probs;
					
					probsum = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
					exprobsum = combine(groupby(probsdf,:household_number), :exprobX => sum)[!,:exprobX_sum];
					
					Xprobk = probsum[household_numbers]; 
					Xexprobk = exprobsum[household_numbers]; 
					 
					gradient_beta[k] =  sum(((1/lambda) .* Xupdated[:,k] .+ ((lambda-1)/lambda) .* Xexprobk .* I_ex_nest .- Xprobk) .* choices .* data_weights .* probsn ./ probabilities);
				end
				
				# Gradient wrt lambda (This should be a scalar!)
				probsdf[!,:probU] = U .* probsn;
				probsdf[!,:exprobU] = U .* ex_probs;
				probUsum = combine(groupby(probsdf,:household_number), :probU => sum)[!,:probU_sum];
				exprobUsum = combine(groupby(probsdf,:household_number), :exprobU => sum)[!,:exprobU_sum];
				Uprob = probUsum[household_numbers];
				Uexprob = exprobUsum[household_numbers];		
				
				probsdf = DataFrame(household_number = household_numbers[uninsured_plans .== 0],
								expvarall = exp.(U[uninsured_plans .== 0]/lambda));
				N_exchange_household = combine(groupby(probsdf,:household_number), :expvarall => sum)[!,:expvarall_sum];
				N_exchange = N_exchange_household[household_numbers];
				
				gradient_lambda = sum((-(1/lambda^2) .* U + (1 .- uninsured_plans) .* log.(N_exchange) +  
								((1 - lambda)/lambda^2) .* Uexprob .* I_ex_nest .-
								inside_option_probs .* log.(N_exchange) + (1/lambda) .* Uprob) .*
								choices .* data_weights .* probsn ./ probabilities);

				gradient += 1/ns .* cat(dims=1,gradient_beta,gradient_lambda);
			end
			
			grad[1:(length(all_covariates) + 1)] = gradient;
		end
	
		function calculate_rp_seq_logit_gradient_for_test(min_beta)
			
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			prob_matrix = zeros(length(choices),ns);
			sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				hseq_probs_updated[!,:prob] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
					household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
				prob_matrix[:,n] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
			end
			sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
		
			Xupdated = copy(X);
			probsdf = DataFrame(household_number = household_numbers);
			gradient = 0.0;
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				
				probsn = vec(prob_matrix[:,n]);
				sequence_probsn = vec(sequence_prob_matrix[:,n]);
				
				probsdf[!,:probX] = Xupdated[:,1] .* probsn;
				hseq_probs_updated[!,:X] = Xupdated[observed_choice_indices,1]; 
				hseq_probs_updated[!,:Xseq_prob] = combine(groupby(probsdf,:household_number), :probX => sum)[!,:probX_sum];
				Xseq1 = combine(groupby(hseq_probs_updated,:household_id), :X => sum)[!,:X_sum];
				Xseq_prob1 = combine(groupby(hseq_probs_updated,:household_id), :Xseq_prob => sum)[!,:Xseq_prob_sum];
				
				gradient += 1/ns .* (transpose(Xseq1 .- Xseq_prob1) * (hseq_weights .* sequence_probsn ./ sequence_probs));
			end
		
			return(gradient)
		end
	
		function calculate_logit_inattention_log_likelihood_for_hessian(min_beta)
			Xupdated = copy(X);
			if random_parameters
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,1]; 
			end
			
			if nested
				probabilities = simulate_nested_probabilities(cat(dims=1,min_beta[1:(premium_shock_index-1)],last(min_beta)),X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			else
				probabilities = simulate_logit_probabilities(min_beta[1:(premium_shock_index-1)],X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			search_probs = calculate_search_probabilities(min_beta[premium_shock_index:(length(min_beta))],Xsearch);
			return((transpose(choices .* entrant_weights) * log.(probabilities))[1] .+ 
					(transpose(choices .* switcher_weights) * log.(search_probs .* probabilities))[1] .+ 
					(transpose(choices .* nonswitcher_weights) * log.((1 .- search_probs) .+ search_probs .* probabilities))[1])
		end
	
		function calculate_logit_log_likelihood_for_hessian(min_beta)
			Xupdated = copy(X);
			prob_matrix = zeros(length(choices),ns);
			for n in 1:ns
				if control_function & use_eta
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
				end
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		
			return((transpose(choices .* data_weights) * log.(probabilities))[1])
		end
	
		function calculate_nested_sequence_log_likelihood_for_test(min_beta)
		
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
			for n in 1:ns
				Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
				hseq_probs_updated[!,:prob] = simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],
					household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
				sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
			end
			sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
		
			return(transpose(hseq_weights) * log.(sequence_probs))
		end
	
	#### Form Random Parameter Covariate Matrix
	
		function form_rp_mean_vector(random_variables)
			mean_vector = [betap[findall(random_variables .== "Anthem")[1]],
							betap[findall(random_variables .== "Blue_Shield")[1]],
							betap[findall(random_variables .== "Kaiser")[1]],
							betap[findall(random_variables .== "Health_Net")[1]],
							betap[findall(random_variables .== "AV")[1]]];
			return(mean_vector)
		end
	
		function form_rp_cov_matrix(betap)
			cov_matrix = [betap[findall(all_covariates .== "Anthem_random")[1]] betap[findall(all_covariates .== "BS_Anthem_cov")[1]] betap[findall(all_covariates .== "Kaiser_Anthem_cov")[1]] betap[findall(all_covariates .== "HN_Anthem_cov")[1]] betap[findall(all_covariates .== "AV_Anthem_cov")[1]];
							betap[findall(all_covariates .== "BS_Anthem_cov")[1]] betap[findall(all_covariates .== "Blue_Shield_random")[1]] betap[findall(all_covariates .== "Kaiser_BS_cov")[1]] betap[findall(all_covariates .== "HN_BS_cov")[1]] betap[findall(all_covariates .== "AV_BS_cov")[1]];
							betap[findall(all_covariates .== "Kaiser_Anthem_cov")[1]] betap[findall(all_covariates .== "Kaiser_BS_cov")[1]] betap[findall(all_covariates .== "Kaiser_random")[1]] betap[findall(all_covariates .== "HN_Kaiser_cov")[1]] betap[findall(all_covariates .== "AV_Kaiser_cov")[1]];
							betap[findall(all_covariates .== "HN_Anthem_cov")[1]] betap[findall(all_covariates .== "HN_BS_cov")[1]] betap[findall(all_covariates .== "HN_Kaiser_cov")[1]] betap[findall(all_covariates .== "Health_Net_random")[1]] betap[findall(all_covariates .== "AV_HN_cov")[1]];
							betap[findall(all_covariates .== "AV_Anthem_cov")[1]] betap[findall(all_covariates .== "AV_BS_cov")[1]] betap[findall(all_covariates .== "AV_Kaiser_cov")[1]] betap[findall(all_covariates .== "AV_HN_cov")[1]] betap[findall(all_covariates .== "AV_random")[1]]];
			return(cov_matrix)
		end
	
		function form_lower_triangular_matrix(betap)
			ll_matrix = [betap[findall(all_covariates .== "Anthem_random")[1]] 0 0 0 0;
							betap[findall(all_covariates .== "BS_Anthem_cov")[1]] betap[findall(all_covariates .== "Blue_Shield_random")[1]] 0 0 0;
							betap[findall(all_covariates .== "Kaiser_Anthem_cov")[1]] betap[findall(all_covariates .== "Kaiser_BS_cov")[1]] betap[findall(all_covariates .== "Kaiser_random")[1]] 0 0;
							betap[findall(all_covariates .== "HN_Anthem_cov")[1]] betap[findall(all_covariates .== "HN_BS_cov")[1]] betap[findall(all_covariates .== "HN_Kaiser_cov")[1]] betap[findall(all_covariates .== "Health_Net_random")[1]] 0;
							betap[findall(all_covariates .== "AV_Anthem_cov")[1]] betap[findall(all_covariates .== "AV_BS_cov")[1]] betap[findall(all_covariates .== "AV_Kaiser_cov")[1]] betap[findall(all_covariates .== "AV_HN_cov")[1]] betap[findall(all_covariates .== "AV_random")[1]]];
			return(ll_matrix)
		end
	
		function transform_std_multivariate_normal(V,ll_matrix)
			#random_variates = transpose(cholesky(cov_matrix).L * transpose(V));
			random_variates = transpose(ll_matrix * transpose(V));
			random_variates[:,findall(random_variables .== "AV")[1]] = exp.(random_variates[:,findall(random_variables .== "AV")[1]]);
			return(random_variates)
		end
	
	### Calculate Objective
	
		function newton_nested_logit(betap,X,choices,tolerance,control_function)
		
			Xupdated = copy(X);
		
			# Initialization
			if random_parameters & sequence_flag
				hseq_probs_updated = copy(hseq_probs);
				prob_matrix = zeros(length(choices),ns);
				sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hseq_probs_updated[!,:prob] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
					prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
					sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
				end
				sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
			elseif random_parameters & control_function
				prob_matrix = zeros(length(choices),ns);
				for n in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
			elseif random_parameters 
				prob_matrix = zeros(length(choices),ns);
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
			else
				probabilities = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			
			obj = calculate_log_likelihood(choices,probabilities,data_weights);
			obj_vector = copy(obj);
			objprev = 0;
			println(obj)
			
			# Now iterate
			while abs(obj - objprev) > tolerance
			
				# Save old objective
				objprev = obj;
				
				# Calculate gradient and hessian
				grad = compute_nested_logit_gradient(betap,Xupdated,choices,probabilities,household_numbers,data_weights);
				hess = compute_nested_logit_hessian(betap,Xupdated,choices,probabilities,household_numbers,data_weights);
				#grad = compute_nested_logit_gradient(betap,Xupdated,choices,probabilities,household_numbers,ones(length(data_weights)));
				#hess = compute_nested_logit_hessian(betap,Xupdated,choices,probabilities,household_numbers,ones(length(data_weights)));
					
				# Update value of beta
				betap = betap - inv(hess) * grad;
				if random_parameters & sequence_flag
					hseq_probs_updated = copy(hseq_probs);
					prob_matrix = zeros(length(choices),ns);
					sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
					for n in 1:ns
						Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
						hseq_probs_updated[!,:prob] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
							household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
						prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
						sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
					end
					sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
				elseif random_parameters & control_function
					prob_matrix = zeros(length(choices),ns);
					for n in 1:ns
						if use_eta
							Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
						end
						Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
						prob_matrix[:,n] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
					end
					probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
				elseif random_parameters
					probabilities = zeros(length(choices));
					for n in 1:ns
						Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
						probabilities += simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
					end
				else
					probabilities = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				obj = calculate_log_likelihood(choices,probabilities,data_weights); 
			
				println(obj)
				writedlm(filename,betap)
			end
			
			return(betap)
		end
	
		function newton_logit(betap,X,choices,tolerance,control_function)
		
			Xupdated = copy(X);
			Vupdated = copy(V);
		
			# Initialization
			if random_parameters & sequence_flag
				hseq_probs_updated = copy(hseq_probs);
				prob_matrix = zeros(length(choices),ns);
				sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
				for n in 1:ns
					if control_function & use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					
					# Transform Standard to Multivariate Normal
					if add_random_cov_terms
						cov_matrix = form_rp_cov_matrix(betap);
						Vupdated[:,:,n] = transform_std_multivariate_normal(V[:,:,n],cov_matrix);
					end
					
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* Vupdated[:,:,n];
					
					hseq_probs_updated[!,:prob] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
					prob_matrix[:,n] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
					sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod];
				end
				sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
			elseif random_parameters
				prob_matrix = zeros(length(choices),ns);
				for n in 1:ns
					if control_function & use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					prob_matrix[:,n] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
			else
				probabilities = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			
			if inattention
				obj = calculate_inattention_log_likelihood(choices,probabilities,search_probs);
			elseif random_parameters & sequence_flag
				obj = calculate_sequence_log_likelihood(sequence_probs,hseq_weights);
			else	
				obj = calculate_log_likelihood(choices,probabilities,data_weights);
				#obj = calculate_log_likelihood(choices,probabilities,ones(length(data_weights)));
			end
			
			objprev = -1e16;
			println(obj)
			
			# Now iterate
			while (abs(obj - objprev) > tolerance) #& (obj > objprev)
			
				# Save old objective
				objprev = copy(obj);
				
				# Update value of beta
				if random_parameters & sequence_flag
					grad = compute_rp_logit_seq_gradient(X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights);
					hess = compute_rp_logit_seq_hessian(X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights);
				elseif random_parameters
					grad = compute_rp_logit_gradient(X,choices,probabilities,prob_matrix,data_weights);
					hess = compute_rp_logit_hessian(X,choices,probabilities,prob_matrix,data_weights);
				else
					grad = compute_logit_gradient(Xupdated,choices,probabilities,data_weights);
					hess = compute_logit_hessian(betap,Xupdated,choices,probabilities,household_numbers,data_weights);
					#grad = compute_logit_gradient(Xupdated,choices,probabilities,ones(length(data_weights)));
					#hess = compute_logit_hessian(betap,Xupdated,choices,probabilities,household_numbers,ones(length(data_weights)));
				end
			
				betap -= inv(hess) * grad;
			
				# Calculate new objective value
				if random_parameters & sequence_flag
					hseq_probs_updated = copy(hseq_probs);
					prob_matrix = zeros(length(choices),ns);
					sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
					for n in 1:ns
						if control_function & use_eta
							Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
						end
						
						# Transform Standard to Multivariate Normal
						if add_random_cov_terms
							cov_matrix = form_rp_cov_matrix(betap);
							Vupdated[:,:,n] = transform_std_multivariate_normal(V[:,:,n],cov_matrix);
						end
						
						Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* Vupdated[:,:,n];
						
						
						hseq_probs_updated[!,:prob] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
							household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
						prob_matrix[:,n] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
						sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod];
					end
					sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
				elseif random_parameters
					prob_matrix = zeros(length(choices),ns);
					for n in 1:ns
						if control_function & use_eta
							Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
						end
						Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
						prob_matrix[:,n] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
					end
					probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
				else
					probabilities = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				
				if inattention
					obj = calculate_inattention_log_likelihood(choices,probabilities,search_probs);
				elseif random_parameters & sequence_flag
					obj = calculate_sequence_log_likelihood(sequence_probs,hseq_weights);
				else	
					obj = calculate_log_likelihood(choices,probabilities,data_weights);
					#obj = calculate_log_likelihood(choices,probabilities,ones(length(data_weights)));
				end
				
				println(obj)
				writedlm(filename,betap)
			end
			
			return(betap)
		end
	
		function objective_nested_logit(betap,X,choices,control_function,random_parameters) 
			writedlm(filename,betap)
			Xupdated = copy(X);
			hseq_probs_updated = copy(hseq_probs);
			
			if control_function & random_parameters & sequence_flag
				sequence_probs = zeros(length(unique(observed_choice_household_ids)));
				for n in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hseq_probs_updated[!,:prob] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices]; 
					sequence_probs += combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
				end
			elseif control_function & random_parameters
				probs = zeros(length(choices));
				for n in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probs += simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif control_function
				probs = zeros(length(choices));
				for m in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,m] ;
					end
					probs += simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif random_parameters & sequence_flag
				sequence_probs = zeros(length(unique(observed_choice_household_ids)));
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hseq_probs_updated[!,:prob] = simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices]; 
					sequence_probs += combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
				end
			elseif random_parameters
				probs = zeros(length(choices));
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probs += simulate_nested_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif inattention
				probs = simulate_nested_probabilities(cat(dims=1,betap[1:(premium_shock_index-1)],last(betap)),X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				search_probs = calculate_search_probabilities(betap[premium_shock_index:(length(betap)-1)],Xsearch);
			else
				probs = simulate_nested_probabilities(betap,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			
			if inattention
				return(-calculate_inattention_log_likelihood(choices,probs,search_probs))
			elseif random_parameters & sequence_flag
				return(-calculate_sequence_log_likelihood(sequence_probs,hseq_weights))
			else	
				return(-calculate_log_likelihood(choices,probs,data_weights))
			end
		end
	
		function objective_logit(betap,X,choices,control_function,random_parameters) 
			writedlm(filename,betap)
			Xupdated = copy(X);
			Vupdated = copy(V);
			hseq_probs_updated = copy(hseq_probs);
			
			if inattention & random_parameters
				probs = zeros(length(choices));
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probs += simulate_logit_probabilities(betap[1:(premium_shock_index-1)],Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
				search_probs = calculate_search_probabilities(betap[premium_shock_index:(length(betap))],Xsearch);
			elseif control_function & random_parameters & sequence_flag
				sequence_probs = zeros(length(unique(observed_choice_household_ids)));
				for n in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					
					# Transform Standard to Multivariate Normal
					if add_random_cov_terms
						ll_matrix = form_lower_triangular_matrix(betap);
						Vupdated[:,:,n] = transform_std_multivariate_normal(V[:,:,n],ll_matrix);
					end
					
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* Vupdated[:,:,n]; 
					
					if smooth_flag
						hseq_probs_updated[!,:prob] = simulate_logit_probabilities_smooth(betap,Xupdated[uninsured_plans .== 0,:],
							household_numbers[uninsured_plans .== 0],uninsured_plans,tau)[observed_choice_indices]; 
					else
						hseq_probs_updated[!,:prob] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
							household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices]; 
					end
					#hseq_probs_updated[!,:prob] = (simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
					#	household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices]) .^ (data_weights[observed_choice_indices]); 
					sequence_probs += combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
				end
			elseif control_function & random_parameters
				probs = zeros(length(choices));
				for n in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probs += simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif control_function
				probs = zeros(length(choices));
				for m in 1:ns
					if use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,m];
					end
					probs += simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif random_parameters & sequence_flag
				sequence_probs = zeros(length(unique(observed_choice_household_ids)));
				for n in 1:ns
					# Transform Standard to Multivariate Normal
					if add_random_cov_terms
						ll_matrix = form_lower_triangular_matrix(betap);
						Vupdated[:,:,n] = transform_std_multivariate_normal(V[:,:,n],ll_matrix);
					end
					
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* Vupdated[:,:,n]; 
					hseq_probs_updated[!,:prob] = simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
					sequence_probs += combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod] ./ ns;
				end
			elseif random_parameters
				probs = zeros(length(choices));
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probs += simulate_logit_probabilities(betap,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif inattention
				probs = simulate_logit_probabilities(betap[1:(premium_shock_index-1)],X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				search_probs = calculate_search_probabilities(betap[premium_shock_index:(length(betap))],Xsearch);
			else
				probs = simulate_logit_probabilities(betap,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
			
			if inattention
				return(-calculate_inattention_log_likelihood(choices,probs,search_probs))
			elseif random_parameters & sequence_flag
				return(-calculate_sequence_log_likelihood(sequence_probs,hseq_weights))
				#return(-calculate_sequence_log_likelihood(sequence_probs,ones(length(hseq_weights))))
			else	
				return(-calculate_log_likelihood(choices,probs,data_weights))
				#return(-calculate_log_likelihood(choices,probs,ones(length(data_weights))))
			end
		end
	
	#### Save choice output
	
	#function save_choice_output(probabilities)
		
		#household_numbers_ins = household_numbers[uninsured_plans .== 0];
	
		# rsid: insurer, HMO, metal, rating area
		
		#choice_output = DataFrame();
		
		#for plan in plans[!,:Plan_Name]
		#	plan_indices = findall(data[!,:plan_name] .== plan);
		#	for covariate in setdiff(product_covariates,["Premium","AV","Anthem_HMO"])
		#		X[plan_indices,findall(all_covariates .== covariate)[1]] .= 
		#			plans[findall(plans[!,:Plan_Name] .== plan)[1],findall(names(plans) .== covariate)[1]];
		#	end
		#end
		
		#choice_output[!,:rating_area] = households[household_numbers_ins,:rating_area];
		#choice_output[!,:metal] .= "Silver";
		#choice_output[data[uninsured_plans .== 0,:av] .== 0.6,:metal] .= "Bronze";
		#choice_output[data[uninsured_plans .== 0,:av] .== 0.8,:metal] .= "Gold";
		#choice_output[data[uninsured_plans .== 0,:av] .== 0.9,:metal] .= "Platinum";
		#choice_output[!,:prob] = probabilities[uninsured_plans .== 0];
		
		#return(choice_output)
	#end
	
	# Simulation exercise to check random parameters estimation is working
	function check_random_parameters(beta_initial,true_solution)
	
		# Get Simulated Choice Probabilities
		Xupdated = copy(X);
		Xupdated[:,findall(all_covariates .== "intercept_random")[1]] = Xupdated[:,findall(all_covariates .== "intercept")[1]] .* V[:,findall(random_product_covariates .== "intercept_random")[1],1];
		Xupdated[:,findall(all_covariates .== "AV_random")] = Xupdated[:,findall(all_covariates .== "AV_random")[1]] .* V[:,findall(random_product_covariates .== "AV_random")[1],1];
		if add_silver
			Xupdated[:,findall(all_covariates .== "Silver_random")] = Xupdated[:,findall(all_covariates .== "Silver_random")[1]] .* V[:,findall(random_product_covariates .== "Silver_random")[1],1];
		end
		Xupdated[:,findall(all_covariates .== "Anthem_random")] = Xupdated[:,findall(all_covariates .== "Anthem_random")[1]] .* V[:,findall(random_product_covariates .== "Anthem_random")[1],1];
		Xupdated[:,findall(all_covariates .== "Blue_Shield_random")] = Xupdated[:,findall(all_covariates .== "Blue_Shield_random")[1]] .* V[:,findall(random_product_covariates .== "Blue_Shield_random")[1],1];
		Xupdated[:,findall(all_covariates .== "Health_Net_random")] = Xupdated[:,findall(all_covariates .== "Health_Net_random")[1]] .* V[:,findall(random_product_covariates .== "Health_Net_random")[1],1];
		Xupdated[:,findall(all_covariates .== "Kaiser_random")] = Xupdated[:,findall(all_covariates .== "Kaiser_random")[1]] .* V[:,findall(random_product_covariates .== "Kaiser_random")[1],1];
		choice_probs_for_sim = simulate_logit_probabilities(true_solution,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
		
		# Get Simulated Choices
		Random.seed!(2);
		
		# Simulate Sequence of Choices
		simulated_plan_choices = households[!,:channel];
		for i in 1:length(unique(observed_choice_household_ids))
			year_in_sequence = 1;
			households_in_sequence = findall(hseq_probs[!,:household_id] .== unique(observed_choice_household_ids)[i]);
			households_in_sequence_years = households[households_in_sequence,:year];
			households_in_sequence = households_in_sequence[sortperm(households_in_sequence_years)];
			
			ids = findall(household_numbers .== households_in_sequence[year_in_sequence]);
			simulated_plan_choices[households_in_sequence[year_in_sequence]] = sample(data[ids,:plan_name],Weights(choice_probs_for_sim[ids]));
			
			while year_in_sequence < length(households_in_sequence)
				year_in_sequence += 1;
				ids = findall(household_numbers .== households_in_sequence[year_in_sequence]);
				
				# Update previous choice variables
				new_prev_choice_index = intersect(findall(household_numbers .== households_in_sequence[year_in_sequence]),
					findall(data[!,:plan_name] .== simulated_plan_choices[households_in_sequence[year_in_sequence-1]]));
				
				if length(new_prev_choice_index) > 1
					new_prev_choice_index = new_prev_choice_index[1];
				
					prev_ind = findall(all_covariates .== "previous_choice")[1];
					
					Xupdated[ids,prev_ind] .= 0;
					Xupdated[new_prev_choice_index,prev_ind] = 1;
					
					Xupdated[ids,findall(all_covariates .== "previous_choice_250to400")[1]] = Xupdated[ids,prev_ind] .* ((households[households_in_sequence[year_in_sequence],:FPL] .> 2.50) & (households[households_in_sequence[year_in_sequence],:FPL] .<= 4)) ; 
					Xupdated[ids,findall(all_covariates .== "previous_choice_gt400")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:FPL] .> 4) ; 
					Xupdated[ids,findall(all_covariates .== "previous_choice_0to17")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_0to17]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_18to34")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_18to25] + households[households_in_sequence[year_in_sequence],:perc_26to34]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_35to54")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_35to44] + households[households_in_sequence[year_in_sequence],:perc_45to54]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_male")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_male]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_family")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:household_size] .> 1);
					Xupdated[ids,findall(all_covariates .== "previous_choice_Asian")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_asian]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_Black")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_black]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_Hispanic")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_hispanic]);
					Xupdated[ids,findall(all_covariates .== "previous_choice_Other")[1]] = Xupdated[ids,prev_ind] .* (households[households_in_sequence[year_in_sequence],:perc_other]);
				
					Xupdated[ids,findall(all_covariates .== "previous_choice_Anthem")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "Anthem")[1]];
					Xupdated[ids,findall(all_covariates .== "previous_choice_Blue_Shield")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "Blue_Shield")[1]];
					Xupdated[ids,findall(all_covariates .== "previous_choice_Kaiser")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "Kaiser")[1]];
					Xupdated[ids,findall(all_covariates .== "previous_choice_Health_Net")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "Health_Net")[1]];			
					Xupdated[ids,findall(all_covariates .== "previous_choice_HMO")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "HMO")[1]];
					Xupdated[ids,findall(all_covariates .== "previous_choice_AV")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "AV")[1]];
					if add_silver
						Xupdated[ids,findall(all_covariates .== "previous_choice_silver")[1]] = Xupdated[ids,prev_ind] .* Xupdated[ids,findall(all_covariates .== "Silver")[1]];
					end
				end
				
				# Recalculate probabilities
			
				choice_probs_for_sim = simulate_logit_probabilities(beta_initial,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				simulated_plan_choices[households_in_sequence[year_in_sequence]] = sample(data[ids,:plan_name],Weights(choice_probs_for_sim[ids]));
			end
		end
		
		sequence_probs = zeros(length(unique(observed_choice_household_ids)));
		
		# Simulate new previous choices
		#simulated_plan_choices = households[!,:channel];
		#for i in 1:size(households)[1]
		#	ids = findall(household_numbers .== i);
		#	simulated_plan_choices[i] = sample(data[ids,:plan_name],Weights(choice_probs_for_sim[ids]));
		#end
		simulated_choices = zeros(size(Xupdated)[1]);
		simulated_choices[findall(in(string.(simulated_plan_choices,"_",collect(1:1:size(households)[1]))),
			string.(data[!,:plan_name],"_",household_numbers))] .= 1;
			
		
		# Resestimate model
		beta_initial_for_sim = copy(beta_initial);
		obj_logit(betap) = objective_logit(betap,X,simulated_choices,control_function,random_parameters);
		logit_results = optimize(obj_logit,convert(Array,beta_initial),iterations=50000,show_trace=true) ;
		
		return(logit_results)
	end
	
	function get_choice_probabilities(min_beta)
		
		Xupdated = copy(X);
		if nested
			if control_function & random_parameters
				probabilities = zeros(length(choices));
				for n in 1:ns
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n] ;
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probabilities += simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif control_function
				probabilities = zeros(length(choices));
				for m in 1:ns
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,m] ;
					probabilities += simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif random_parameters
				probabilities = zeros(length(choices));
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					probabilities += simulate_nested_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
			elseif inattention 
				probabilities = simulate_nested_probabilities(cat(dims=1,min_beta[1:(premium_shock_index-1)],last(min_beta)),X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			else
				probabilities = simulate_nested_probabilities(min_beta,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
			end
		else
			if inattention
				prob_matrix = zeros(length(choices),ns);
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					prob_matrix[:,n] = simulate_logit_probabilities(min_beta[1:(premium_shock_index-1)],Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans) ./ ns;
				end
				return(prob_matrix)
			elseif random_parameters & sequence_flag
				hseq_probs_updated = copy(hseq_probs);
				prob_matrix = zeros(length(choices),ns);
				sequence_prob_matrix = zeros(length(unique(observed_choice_household_ids)),ns);
				for n in 1:ns
					if control_function & use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hseq_probs_updated[!,:prob] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],
						household_numbers[uninsured_plans .== 0],uninsured_plans)[observed_choice_indices];
					prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
					sequence_prob_matrix[:,n] = combine(groupby(hseq_probs_updated,:household_id), :prob => prod)[!,:prob_prod];
				end
				return(cat(dims = 1,prob_matrix,sequence_prob_matrix))
			elseif random_parameters
				prob_matrix = zeros(length(choices),ns);
				for n in 1:ns
					if control_function & use_eta
						Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
					end
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				end
				return(prob_matrix)
			else
				probabilities = simulate_logit_probabilities(min_beta,X[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
				return(probabilities)
			end
		end
		
		return(probabilities)
	end
	
	function calculate_instrumented_probs(min_beta)
		Xupdated = copy(X);
		Xupdated[:,findall(all_covariates .== "Premium")[1]] .-= Xupdated[:,findall(all_covariates .== "residual")[1]];
		prob_matrix = zeros(length(choices),ns);
		for n in 1:ns
			if control_function & use_eta
				Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n];
			end
			Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
			prob_matrix[:,n] = simulate_logit_probabilities(min_beta,Xupdated[uninsured_plans .== 0,:],household_numbers[uninsured_plans .== 0],uninsured_plans);
		end
		probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		return(probabilities)
	end
	
	function process_parameters(min_beta,probabilities);
		
		# Get standard errors
		Xupdated = copy(X);
		if nested
			if control_function & random_parameters
				hess = zeros(length(all_covariates)+1,length(all_covariates)+1);
				for n in 1:ns
					Xupdated[:,findall(all_covariates .== "eta")[1]] = eta[:,n] ;
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hess += compute_nested_logit_hessian(min_beta,Xupdated,choices,probabilities,household_numbers,data_weights) ./ ns;
				end
				se_beta = sqrt.(diag(inv(-hess)));
			elseif random_parameters
				hess = zeros(length(all_covariates)+1,length(all_covariates)+1);
				for n in 1:ns
					Xupdated[:,random_parameter_indices] = X[:,random_parameter_indices] .* V[:,:,n]; 
					hess += compute_nested_logit_hessian(min_beta,Xupdated,choices,probabilities,household_numbers,data_weights) ./ ns;
				end
				se_beta = sqrt.(diag(inv(-hess)));
			else
				se_beta = sqrt.(diag(inv(-compute_nested_logit_hessian(min_beta,Xupdated,choices,probabilities,household_numbers,data_weights))));
			end
		else
			if inattention
				# NOTE: running the Hessian function takes about 8-9 hours!
				hess = hessian(calculate_logit_inattention_log_likelihood_for_hessian,min_beta);
				se_beta = sqrt.(max.(diag(inv(-hess)),0));
			elseif random_parameters & sequence_flag
				se_beta = sqrt.(max.(diag(inv(-compute_rp_logit_seq_hessian(X,choices,sequence_probs,sequence_prob_matrix,prob_matrix,hseq_weights))),0));
			elseif random_parameters
				se_beta = sqrt.(max.(diag(inv(-compute_rp_logit_hessian(X,choices,probabilities,prob_matrix,data_weights))),0));
			else
				se_beta = sqrt.(max.(diag(inv(-compute_logit_hessian(min_beta,Xupdated,choices,probabilities,household_numbers,data_weights))),0));
			end
		end
	
		# Save output
		param_output = DataFrame();
		if nested
			if inattention
				param_output[!,:Variable] = cat(dims=1,all_covariates,"lambda",search_covariates);
			else
				param_output[!,:Variable] = cat(dims=1,all_covariates,"lambda");
			end
		else
			if inattention
				param_output[!,:Variable] = cat(dims=1,all_covariates,search_covariates);
			else
				param_output[!,:Variable] = all_covariates;
			end
		end
		param_output[!,:Coef] = min_beta;
		if random_parameters
			#param_output[random_parameter_indices,:Coef] = abs.(param_output[random_parameter_indices,:Coef]);
		end
		param_output[!,:SE] = se_beta;
		param_output[!,:pval] = 2 * (1 .- cdf.(Normal(),abs.(min_beta ./ se_beta)));
		
		return(param_output)
	end
	
	function compute_elasticities(min_beta,probabilities)

		if (rich_flag | poor_flag | old_flag | young_flag)
			fields = ["Overall","Female","Male","Asian","Black","Hispanic","Other","White","Single","Family"];
		elseif add_complex_interactions
			fields = ["Overall","0-250_0to34","250-400_0to34","400+_0to34","0-250_35to54","250-400_35to54","400+_35to54","0-250_55+","250-400_55+","400+_55+",
				"Female","Male","Asian","Black","Hispanic","Other","White","Single","Family"]; 
		else
			fields = ["Overall","0-250","250-400","400+","Female","Male","0-17","18-34","35-54","55+","Asian","Black","Hispanic","Other","White","Single","Family"];
		end

		elasticity_output = create_elasticity_table(min_beta,X,fields); 
		
		return(elasticity_output)
	end	
	
	function compute_switching_costs(min_beta,probabilities)
		
		if (rich_flag | poor_flag | old_flag | young_flag)
			fields = ["Overall","Female","Male",
								"Asian","Black","Hispanic","Other","White",
								"Single","Family",
								"Anthem","Blue_Shield","Kaiser","Health_Net","Other_Insurer",
								"Non-HMO","HMO"];
		elseif add_complex_interactions
			fields = ["Overall","0-250_0to34","250-400_0to34","400+_0to34","0-250_35to54","250-400_35to54","400+_35to54","0-250_55+","250-400_55+","400+_55+",
				"Female","Male","Asian","Black","Hispanic","Other","White","Single","Family"];
		else
			fields = ["Overall","0-250","250-400","400+","Female","Male",
								"0-17","18-34","35-54","55+",
								"Asian","Black","Hispanic","Other","White",
								"Single","Family",
								"Anthem","Blue_Shield","Kaiser","Health_Net","Other_Insurer",
								"Non-HMO","HMO"];
		end
			
		inertia_output = create_inertia_table(min_beta,X,fields) 
		
		return(inertia_output)
	end
	
	function compute_average_params(min_beta,probabilities)
		
		# Calculate average premium and inertia parameters
		alpha = calculate_alpha(min_beta,X,"All"); 
		if !inattention
			inertia_param = calculate_inertia_param(min_beta,X,"All"); 			
		end
		
		avg_param_output = DataFrame();
		if inattention
			avg_param_output[!,:Parameter] = ["Premium"];
			avg_param_output[!,:Average] = [sum(alpha .* probabilities)./sum(probabilities)];
		else
			avg_param_output[!,:Parameter] = ["Premium","Previous_Choice"];
			avg_param_output[!,:Average] = cat(dims=1,sum(alpha .* probabilities)./sum(probabilities),sum(inertia_param .* probabilities)./sum(probabilities));
		end
		
		return(avg_param_output)
	end